#A lexander F. Gazmararian
# afg2@princeton.edu

# Load packages
library(tidyverse)
library(margins)
library(modelsummary)
library(emmeans)
library(sandwich)
library(estimatr)
# Packages for structural topic model
library(stm)
library(geometry)
library(Rtsne)
library(rsvd)

# Load functions
source("code/fun/book_theme.r")
source("code/fun/savefig.r")
source("code/fun/fix_txt.r")

# Set up coefficient names
source("code/fun/coefnames4tables.r")
coefnames <-
  c(
    coefnames,
    "lock_second" = "Lock-in order",
    "delegate_lead_treat"="Community delegation",
    "delegate_voice_treat"="Public input",
    "delegate_fund_treat"="Long-term funding",
    "delegate_lead_treat:delegate_voice_treat" = "Community delegation x Public input",
    "delegate_lead_treat:delegate_fund_treat" = "Community delegation x Long-term funding",
    "delegate_voice_treat:delegate_fund_treat" = "Public input x Long-term funding",
    "delegate_lead_treat:delegate_voice_treat:delegate_fund_treat" = "Community delegation x Public input x Long-term funding",
    "lockin_treat_poolGeneral Benefits" = "General benefits treatment",
    "lockin_treat_poolSpecific Benefits" = "Specific benefits treatment",
    "consensus_treat" = "Consensus treatment",
    "consensus_treat:PartySummaryNeither" = "Consensus treatment x Independent",
    "consensus_treat:PartySummaryRepublican" = "Consensus treatment x Republican"
  )

# Load data
g <- readRDS("data/NatQual_Winter22.rds")

# Specify baseline model
f.base <- y ~ age + Female + Black + Hispanic +  PartySummary + CollegeDegree + employfull + Rural_bin + Investment + ClimateImpact_tri + income5 + trust_index

# Create function to estimate treatment effects for all outcomes
estimate_delegate <- function (formula) {
  m <- list()
  m[["index"]] <- lm(update(formula, Delegate_index ~ .), g)
  m[["regsupport_bin"]] <- lm(update(formula, RegSupport_bin ~ .), g)
  m[["regsupport_scale"]] <- lm(update(formula, RegSupport_scale ~ .), g)
  m[["benefits_index"]] <- lm(update(formula, DelegateBenefits_index ~ .), g)
  m[["train_bin"]] <- lm(update(formula, Train_bin ~ .), g)
  m[["train_scale"]] <- lm(update(formula, Train_scale ~ .), g)
  m[["localjobs_bin"]] <- lm(update(formula, LocalJobs_bin ~ .), g)
  m[["localjobs_scale"]] <- lm(update(formula, LocalJobs_scale ~ .), g)
  m[["fund_bin"]] <- lm(update(formula, Fund_bin ~ .), g)
  m[["fund_scale"]] <- lm(update(formula, Fund_scale ~ .), g)
  m[["reverse_bin"]] <- lm(update(formula, AssistReverse_bin ~ .), g)
  m[["reverse_scale"]] <- lm(update(formula, AssistReverse_scale ~ .), g)
  m[["input_index"]] <- lm(update(formula, DelegateInput_index ~ .), g)
  m[["trust_bin"]] <- lm(update(formula, DelegateTrust_bin ~ .), g)
  m[["trust_scale"]] <- lm(update(formula, DelegateTrust_scale ~ .), g)
  m[["input_bin"]] <- lm(update(formula, InterpInput_bin ~ .), g)
  m[["input_scale"]] <- lm(update(formula, InterpInput_scale ~ .), g)
  return(m)
}

# Specify model
f.delegate <- update(f.base, y ~ delegate_lead_treat * delegate_voice_treat * delegate_fund_treat + . + delegate_count)

# Estimate models
delegate <- estimate_delegate(f.delegate)

## Create online appendix tables ----

# Extract models for the scale outcome
scale.out <- list(delegate$regsupport_scale, delegate$trust_scale, delegate$input_scale, delegate$reverse_scale, delegate$train_scale, delegate$localjobs_scale)

# Create names for models
names(scale.out) <- c("Support", "Trust", "Input", "Reverse", "Train", "Local")

# Create table output
file <- "tables/ch5/ols_delegation.txt"
modelsummary(
  scale.out,
  vcov = "HC2",
  stars = c("*"=.1,"**"=.05,"***"=.01),
  coef_map = coefnames,
  gof_map = c("nobs", "adj.r.squared"),
  escape=FALSE,
  output="latex",
  ) %>%
  cat(., file = file)
fix_txt(file)

# Extract models for the binary outcome
bin.out <- list(delegate$regsupport_bin, delegate$trust_bin, delegate$input_bin, delegate$reverse_bin, delegate$train_bin, delegate$localjobs_bin)

# Create names for models
names(bin.out) <- c("Support", "Trust", "Input", "Reverse", "Train", "Local")

# Create table output
file <- "tables/ch5/ols_delegation_bin.txt"
modelsummary(
  bin.out,
  vcov = "HC2",
  stars = c("*"=.1,"**"=.05,"***"=.01),
  coef_map = coefnames,
  gof_map = c("nobs", "adj.r.squared"),
  escape=FALSE,
  output="latex",
) %>%
  cat(., file = file)
fix_txt(file)

## Create Figure 5.8 ----

# Function to retrieve average marginal treatment effects
delegate.ame <- lapply(delegate, function(x) {
  ame <-
    margins::margins_summary(
      x,
      variables = c("delegate_lead_treat", "delegate_voice_treat", "delegate_fund_treat"),
      vcov = sandwich::vcovHC(x, "HC2")
    )
  ame <- data.frame(ame)
  return(ame)
})

# Extract outcomes
delegate.df <- bind_rows(delegate.ame, .id = "outcome")

# Construct confidence intervals
delegate.df$lower90 <- with(delegate.df, AME + SE * qnorm(.05))
delegate.df$upper90 <- with(delegate.df, AME + SE * qnorm(.95))

# Function to process the names of outcomes
prep_names <- function(x) {
  x <-
    x %>%
  dplyr::mutate(
    outcome = dplyr::case_when(
      outcome == "regsupport_scale" ~ "Support for regulation",
      outcome == "trust_scale" ~ "Trust in program administration",
      outcome == "input_scale" ~ "Input in process",
      outcome == "reverse_scale" ~ "Durability of assistance"
    ),
    outcome = factor(
      outcome,
      ordered = TRUE,
      levels = c(
        "Durability of assistance",
        "Input in process",
        "Trust in program administration",
        "Support for regulation"
      )
    ),
    factor = dplyr::case_when(
      factor == "delegate_fund_treat" ~ "Long-term funding",
      factor == "delegate_lead_treat" ~ "Delegation to community",
      factor == "delegate_voice_treat" ~ "Public input"
    )
  )
  return(x)
}

# Create plot for book
plot.delegate <-
  delegate.df %>%
  filter(outcome %in% c("regsupport_scale","trust_scale",  "input_scale","reverse_scale")) %>%
  prep_names() %>%
  mutate(
    outcome = factor(
      outcome,
      ordered = TRUE,
      levels = rev(
        c(
          "Durability of assistance",
          "Input in process",
          "Trust in program administration",
          "Support for regulation"
        )
      )
    )
  ) %>%
  ggplot(aes(x = factor, y = AME, color = factor, shape = factor)) +
  geom_hline(yintercept = 0, lty = "dashed", color = "grey") +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0, position = position_dodge(.4)) +
  geom_errorbar(aes(ymin = lower90, ymax = upper90), width = 0, size = 1.75, position = position_dodge(.4)) +
  geom_point(size = 4, position = position_dodge(.4)) +
  scale_x_discrete(label = scales::label_wrap(15)) +
  scale_y_continuous(expand = c(0,0)) +
  facet_wrap(~ outcome) +
  scale_color_grey() +
  book_theme +
  labs(color="Treatment:",shape="Treatment:",x="Treatment",y="Estimate", caption="Notes: Thin and thick lines denote 95 and 90 percent confidence intervals") +
  theme(legend.position="",legend.box.margin = margin(l=-75))
plot.delegate
savefig(plot.delegate, "5.8_figure_delegate", height = 3, filepath = "figures/")

## Examine interaction effect of long-term funding and public input ----

# Result reported around pp. 139-142.
cate.voice <- margins_summary(delegate$input_scale, variables = c("delegate_fund_treat"), at = list(delegate_voice_treat=0:1), vcov = vcovHC(delegate$input_scale, "HC2"))

# Create plot for online appendix
plot.cate <-
  cate.voice %>%
  mutate(
    factor = ifelse(factor == "delegate_fund_treat", "Long-term funding"),
    delegate_voice_treat = ifelse(delegate_voice_treat == 1, "Yes", "No")
  ) %>%
  ggplot() +
  geom_vline(xintercept = 0, color = "grey", lty = "dashed") +
  geom_pointrange(aes(y = factor, x = AME, xmin = lower, xmax = upper, color = delegate_voice_treat, shape = delegate_voice_treat),
                position = position_dodge(.5), fatten = 6) +
  book_theme +
  scale_color_grey() +
  labs(x = "Average Marginal Effect of Long-Term Funding", y = "", color = "Public Input", shape = "Public Input", title = "Outcome: Perceptions of Local Input") +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )
savefig(plot.cate, filename = "si_fig_fundvoice_inter", filepath = "figures/ch5/", height = 2)

# Funding
margins_summary(delegate$input_scale, variables = c("delegate_voice_treat","delegate_lead_treat"), at = list(delegate_fund_treat=0:1), vcov = sandwich::vcovHC(delegate$input_scale, "HC2"))

# Community delegation
margins_summary(delegate$input_scale, variables = c("delegate_voice_treat","delegate_fund_treat"), at = list(delegate_lead_treat=0:1), vcov = vcovHC(delegate$input_scale, "HC2"))

## Structural topic model of open-ended responses ----

# Prepare text
ProcDelgTop <- textProcessor(g$DelegationOE, metadata = g)
DelgTop <- prepDocuments(ProcDelgTop$documents, ProcDelgTop$vocab, ProcDelgTop$meta)

# Estimate model
DelgTop.stm <- stm(
  documents = DelgTop$documents,
  vocab = DelgTop$vocab,
  K = 0,
  seed = 20230115,
  prevalence = ~ delegate_lead_treat + delegate_voice_treat + delegate_fund_treat + delegate_count,
  data = DelgTop$meta,
  init.type = "Spectral"
)

# Create plot of topics
png("figures/ch5/si_delegate_stm.png", width = 6.5, height = 7, units = "in", res = 300, pointsize = 10)
plot(DelgTop.stm, type = "summary")
dev.off()

# Inspect most frequent topics
findThoughts(DelgTop.stm, texts = DelgTop$meta$DelegationOE, n = 45, topic = 35)#mention the government,especially distrust

# Estimate model
DelgTop.effect <- estimateEffect(
  formula = c(35) ~ delegate_lead_treat + delegate_voice_treat + delegate_fund_treat,
  DelgTop.stm,
  meta = DelgTop$meta,
  uncertainty = "Global"
  )

m.deleg <- summary(DelgTop.effect, topics = 35)
m.deleg.df <- data.frame(m.deleg$tables[[1]])
m.deleg.df <- rownames_to_column(m.deleg.df, var = "term")
names(m.deleg.df) <- c("term", "estimate", "std.error", "t", "p.value")
m.deleg.df
# Create N data frame
gl <- data.frame(N = nrow(DelgTop$meta)[[1]])
# Create list
mod <- list(tidy = m.deleg.df, glance = gl)
class(mod) <- "modelsummary_list"
# Create table
file <- "tables/ch5/tbl_stm_delegate.txt"
modelsummary(
  mod,
  stars = c("*"=.1,"**"=.05,"***"=.01),
  coef_map = coefnames,
  output = "latex"
) %>%
  cat(., file = file)
fix_txt(file)
